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Abstract. We show that information on the probability density of local fluctuations 
can be obtained from a numerical renormalisation group calculation of a reduced 
density matrix. We apply this approach to the Anderson-Holstein impurity model 
C^ ' to calculate the ground state probability density p{x) for the displacement x of the 

local oscillator. From this density we can deduce an effective local potential for the 

, oscillator and compare its form with that obtained from a semiclassical approximation 

pH . as a function of the coupling strength. The method is extended to infinite dimensional 

O I Holstein-Hubbard model using dynamical mean field theory. We use this approach to 

I ^i ' compare the probability densities for the displacement of the local oscillator in the 

normal, antiferromagnetic and charge ordered phases. 
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Probability distribution of fluctuations in Anderson- Holstein and Holstein- Hubbard Models2 

1. Introduction 

The numerical renornialization group [H E] (NRG) approacli lias been successfully 
applied to the calculation of static and dynamical response functions for models of 
magnetic impurities and quantum dots, and also been applied to lattice models, such as 
the Hubbard model, in the framework of dynamical mean field theory (DMFT) [3]. The 
response functions calculated give information about the low energy fluctuations which 
are of particular interest in the strong correlation regime. However, as a thermodynamic 
average is taken in calculating these response functions, some of the information about 
these fluctuations, which is contained in the original many-body states, is lost. Here we 
show that if we take only a partial thermodynamic average, such as in the calculation of 
a reduced density matrix, we can learn more about the nature of the local fluctuations, 
and how they vary as a function of the interaction terms. 

We illustrate the approach first of all for the Anderson-Holstein impurity model. 
In the Anderson-Holstein model the occupation of the impurity state is linearly coupled 
to a local harmonic oscillator, which has spatial coordinate x, representing the lattice 
degrees of freedom around the impurity site. As the coupling of the oscillator to the 
impurity state is increased, the nature of the local fluctuations of the oscillator changes. 
If we take a full thermodynamic average then we only get averaged information about 
the oscillator. If we take a partial thermodynamic average, and calculate the reduced 
density matrix at the impurity site, treating the rest of the system as an 'environment', 
we can learn more about how the local lattice fluctuations vary as the coupling strength 
increases. From the reduced density matrix, calculated using the NRG, we can deduce 
the probability distribution p{x) for the x coordinate of the oscillator. From p{x) we 
can also deduce an effective potential Vcs{x), and study the change of this potential 
as a function of the interaction strength and of the frequency of the local oscillator. 
Later in the paper, we extend the method to the infinite dimensional Holstein-Hubbard 
lattice model using dynamical mean field theory (DMFT) in combination with the NRG. 
We then compare how the local probability distribution p{x) changes in the normal, 
antiferromagnetic and charge ordered states for this model. 

2. Local Fluctuations in the Anderson-Holstein Impurity Model 

The Anderson-Holstein model corresponds to the single impurity Anderson model [1] 
with an additional linear coupling to a local phonon mode, as in the Holstein model [5]. 
The Hamiltonian takes the form, 

H = ^ Effif^^ + Uhf^^fif^^ + g{b^ + b){^ fif^^ - 1) 

a (J 

+ Y. M^l.c^a + h.c.) + J2 ^fccLcfe. + (^ob^b. (1) 

k,a k(j 

The impurity level ej, as in the usual Anderson model, is hybridised with conduction 
electrons of the host metal via a matrix element Vfe, with an interaction term U between 
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the electrons in the locahzed f (or d) state. There is in addition a couphng g of 
the impurity site occupation rijo- = CrCr^ to a local oscillator of frequency coq. A 
measure of the hybridisation is the width factor, A(a;) = vr^^ V^6{uj — Ek), which for 
a flat conduction band of width W = 2D, and Vj- independent of fc, we can take as a 
constant A = itV'^/2D. The oscillator coordinate operator x in terms of the creation 
and annihilation operators, b'' and b, is given by a; = (6 + b^)/^/2u^, where we have 
taken the mass of the oscillator, M = 1. We have also set h = 1 so that x^^ has 
the dimension of of square root of energy, and it is convenient to define characteristic 
length scale by Xq = 1/y/uJo- A convenient measure of the effects of the electron-phonon 
coupling on the electronic system is the parameter A = 2g'^ /loq. In the limit Wq — ?■ oo, 
such that A remains finite, the model maps into the Anderson model with f/ — )■ f/ — A, 
Ef ^ Ef + A/2. The behaviour of the model has been studied using the NRG [6j and 
it has been used to study the transport through a quantum dot in the presence of a 
coupling of the occupation dot to local phonon modes [7J . 

To learn more about the state of the oscillator we use the NRG to calculate the 
reduced density matrix prcd at the impurity site. A procedure for calculating the reduced 
density matrix was introduced into NRG calculations by Hofstetter [S]. The original 
motivation was to find an improved way of calculating the higher energy features in 
the spectral density of the impurity Green's function in cases of broken symmetry. For 
NRG calculations the system is recast in the form of a linear chain with impurity at 
one end. Sequential diagonalisations are then carried out starting at the impurity site. 
The information from the shorter length chains is used to calculate the higher energy 
features in the spectral density, and the longer chains the low energy features. In 
the case of broken symmetry the ground state for the shorter chains underestimates the 
degree of symmetry breaking. In Hofstetter's modified procedure the ground state is first 
estimated from the longest chain calculated, and then used to deduce the density matrix 
of the sites corresponding to shorter chain lengths. Incorporating the reduced density 
matrix into the calculation of the spectral density from the shorter chains corrects the 
deficiencies of the standard approach in the case of broken symmetry. Refinements of 
this approach have been introduced more recently based on the use of a complete set of 
NRG states |9l [10]. These have the advantage that the sum rules on the total spectral 
density are satisfied exactly, rather than approximately as in earlier versions of the NRG 
approach. 

The calculation of the density matrix gives additional information which we can 
exploit. For example in the Anderson-Holstein model, if we work backwards from the 
longest chain, we can deduce the reduced density matrix at the impurity site Pf, red- The 
matrix elements of this reduced density matrix will be with respect to the basis states at 
the impurity site, the states of all the other sites are averaged over as they are taken to 
be part of the environment. The electronic states at the impurity site can be labeled by 
the local charge q {q = '^cr^f,o-)y ^^^ the 2;-component of spin m^, and the index of the 
harmonic oscillator states z/. The matrix will be diagonal with respect to the spin and 
charge indices, and so a typical matrix element can be expressed as {pf^rediq, ^tt- z))u,u'- 
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The probability distribution function p{x) for the oscillator coordinate x is given by 

p{x) = '^p{x:q,m^) (2) 



q,m^ 



where 



p{x : g,m^) = y^ 0^(x)(p/,rcd(g, m^))^y(f)l,{x), (3) 



and (j),y{x) is the normalised real space harmonic oscillator wavefunction, 

"^^^""^ = (;^) ^"'^"^'^'^'^(v^^)' (4) 

with the Hermite polynomial Hy{x) of order v. 

It is possible to define an effective potential for the oscillator using an effective 
wavefunction defined by ^eff(2;) = \/p{x), which is taken to be a solution of the one- 
dimensional Schrodinger equation (recall M = 1, h = 1), 

_ IdSpcfiix) 
2 dx^ 
or 



Vcs{x)4Jes{x) = Eipeslx), (5) 



t4.W = ^ + i?#- W 



2^cff(a; 
In terms of p(x) this translates into 



V,six) =E + ^ 



' p"{x) 1 fp\xy^' 



p{x) 2 \ p{x 



(7) 



By construction this potential is such that the ground state wavefunction can be used 
to reproduce the NRG derived p{x). 

2.1. Results for the symmetric model 

Unless otherwise stated we use in the following the parameters W = 2D = 2, vrA = 0.1 
and c^o = 0.1 setting the energy scales for electrons and phonons in this section. The 
phonon frequency cuq has been chosen so that 1/wo is on the scale of the life time of a 
electron on the impurity site, cjq ~ vrA. We know in the adiabatic limit uq ^ ttA, the x 
coordinate becomes a classical variable, and in the opposite limit cjq — > oo it maps on to 
an effective Hubbard model, so the range uq ~ ttA is the interesting one to investigate. 
In Fig. [1] we give results for the probability distribution function p{x) for the 
symmetric model calculated as outlined above. 

The results on the left hand panel correspond to U = {sd = 0). We take as a 
relative measure of the strength of the phonon coupling and hybridisation scale the 
dimensionless factor a = A/ttA ranging from weak {a ^ 1) to strong coupling (a ^ 1). 
As a increases the distribution broadens, and for a = 1.5 a two peak structure can be 
seen which becomes more marked as the coupling strength is increased further. 
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Figure 1. The total probability distribution function p{x) for the oscillator 
displacement x in the ground state for a range of values of a, in the left panel with 
t/ = and the right with [//ttA = 2. 




Figure 2. The components, p{x : 0, 0) (dot-dashed curve), p{x : 1, 1/2)+ p{x : 1, —1/2) 
(full curve) and p{x : 2,0) (dashed) oi p{x) (dots) for the case a — 0.5 (left panel), and 
a = 2.0 (right panel) for U — Q. 



In Fig. [2] we plot the individual components of p(x); p{x : 0, 0), Yl± P{^ '■ !> il/2) 
and p{x : 2,0), corresponding to g = 0,1,2, for U = 0, for the two cases a = 0.5 
(left panel) and a = 2 (right panel). One can see from these curves the two factors 
that lead to the two peak structure on increasing a. One factor is that the maxima 
of curves corresponding to g = and q = 2 are shifted on either side of central peak 
corresponding to g = 1. This reflects that fact that in these charge states for an isolated 
impurity the oscillator is displaced from a; = to vA/^o and —yX/uo respectively. The 
other factor is that the weights of the peaks at g = and g = 2, compared to the weight 
of the central peak corresponding to g = 1, increases with a. The integrated weight 
under the curve Pg is a measure of the probability of the occupation of the local level 
in the state with charge quantum number g. This shift in relative weights is due to the 
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fact that the couphng to the phonon mode induces a local attraction. The weight P2 
measures the probability of the local level being doubly occupied which is equal to the 
expectation value (?^/,t'^'/,4,)• Fot a = 2, the weights Pq = P2 = 0.430261 are in precise 
six figure agreement with (nf^^nf^i) as calculated directly from the NRG calculation, 
and we have Pi = 1 — Pq — P2 = 0.139478. Notice that for a = the corresponding 
values are Pq = P2 = 0.25, Pi = 0.5. The relative shift of the weights can be explained 
to a large extent, but not completely, by the local induced attractive interaction. If we 
take a model with a local attractive interaction U/ttA = —2, but no phonon coupling, 
then the values deduced from the NRG for Pg are Pq = P2 = 0.38289, Pi = 0.23422, 
which underestimates the relative shift away from the state q = 1 found in the phonon 
coupled model with a = 2. This shows that the energetic gain of creating zero and 
doubly occupied sites is higher in the system with phonons as compared to the case 
with instantaneous attraction. 



PfCai) 6 




Figure 3. The spectral density Pf{uj) for tlie impurity single electron Green's function 
for [/ = and a = 1, 1.5. 



The spectral density for the local f-electron Green's function Pf{uj) is shown in Fig. 
|3]for U = ioT a = 1 and a = 1.5. It can be seen that the development of the two peak 
form in p{x) correlates with the development of three peaks in Pf{u}). Clear shoulders 
can be seen in the result for Pf{uj) for a = 1, and the three peak form has fully emerged 
for a = 1.5. The central peak becomes extremely narrow for values of a in the strong 
coupling regime [6j. 

The effect of the interaction U in suppressing the onset of the double peak 
distribution can be seen by comparing the U = results with those shown in the 
right panel of Fig. [1] which are for the case U/nA = 2 [sf = —U/2). The double peak 
structure, which develops in the [7 = case for a ~ 1.5, develops in the U/tiA = 2 
case only when a ~ 2.5. This is a smaller value than what is expected naively in the 
tJo — ^ 00 case where Ues = U — X signaling again that for finite uq when phonons can 
be excited the occupation of the zero and doubly occupied sites is energetically more 
favourable. 

The corresponding effective potentials Ves{x) for the U = case, as deduced from 
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Figure 4. The effective potential Vcs{x) for [/ = and a range of values of a, 
corresponding to the results for p(x) in the left panel of figure [T] 



equation ([7]), are shown in Fig. |H The on-set of a double well feature can be seen 
in the results for a = 1.0, which is before a double peak structure can be seen in the 
corresponding p{x). The effective potentials for higher values of a have a clear double 
well form. 

It is of interest to compare this potential with one calculated using a semiclassical 
approximation, in which we neglect the kinetic energy of the oscillator and treat the 
coordinate a; as a classical variable. This is a commonly used approximation in taking 
the electron-phonon coupling into account, and corresponds to the Born-Oppenheimer 
approximation. Evaluating the impurity contribution to the total ground state energy 
E{x) as a function of x one arrives at an expression for the effective semiclassical 
potential, K-ci(a;) (= E{x)) given by 

2A cu^x^ 2£, 



K--ci(a;) 



'^^%an-'f'^^^^ 



A, 

— log 

It 



e}{x) 



+ A' 



D^ 



ix) with Ve?t{x) deduced from 



TT 2 7r """" V A 

where ef{x) = Sf + yj^uj^gx. 

In Fig. |5]we compare the semiclassical potential Vs. 
equation ([7]) for a = 1 and a = 2. 

The relevant comparison is in the shapes of these potentials, not their absolute 
values, and they have been subject to a constant shift so their forms can be compared 
more easily. It can be seen that the potentials develop in a similar way as the coupling 
strength is increased, though when the double well form develops the potential barrier 
between the wells is less marked in the semiclassical case. The minima occur at very 
similar x/xo-values; —0.74 (V^_ci(x)) and —0.79 (V^ff(a:)) for a = 1, and —1.26 (Vi_ci(x)) 
and —1.35 (Veff(x)) for a = 2. 

The coefficient of the x^ term in Vs_ci(x) changes sign for a = 0.5, which is the 
point at which the double well begins to form. For a > 0.5 there are two mean field 
solutions for the expectation value (x) corresponding to the two minima in Vs-^c\{x) and 
the local symmetry is broken. The positions of the minima in V(_.ci(a^) can be deduced 



Probability distribution of fluctuations in Anderson- Holstein and Holstein- Hubbard ModelsS 



v,„fr) 




■ a=l 

■ a = 2 



-1 1 

x/x„ 



Figure 5. The effective potentials Vcs{x) (full lines) as deduced from equation (O, 
compared with corresponding the semiclassical potentials V^_ci(x) (dashed lines), Eq. 
dHl), for C/ = 0. 



from the equation, 

9v;_ei(x) 



dx 



= UqX + y/2uJog{nf — 1) 



where Uf is the mean field occupation value given by 



. 2 1 fef(x) 



yielding 



X = {x) 



MF 



1 



Uf-l) 



(9) 



(10) 



:ir 



We can deduce an exact relation of this type by introducing an additional to term of the 
form c{b + b^) into the Hamiltonian. The ground state energy E{c) will be a function of 
the coupling c introduced, and we can deduce that 

dE{c) 



dc 



{b + b^) = V2U'o{x). 



(12) 



c=0 



If we now perform a canonical transformation H' = U ^HU with U = e ^^^ ^)l^o ^ the 
terms in H' which depend on c are 



c^ 2^c 

Wo Wo 



(13) 



As the canonical transformation does not effect the energy values we can use this result 
in equation (IT2|) to determine {b + b^), which leads to the result. 



(x) 



uo 



i{nf)-l). 



(14) 



This is the same formula as derived from the semi-classical approximation as given in 
equation (ITT]) , except that the mean field value Uf is replaced by the exact value (nj) 
for the impurity occupation. This result holds for U ^ and is exact. 
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In the mean field broken symmetry solutions we have the two broken symmetry 
solutions at strong coupling, such that rij ~ and Uf ~ 2, giving from equation (ITT]) 
(£)mf ~ ivA/cJo for the positions of the two minima V^_ci(a;). In the exact solution for 
the impurity model, however, this broken symmetry must be restored and the average 
value of X must be zero. This is clearly the case in the NRG solution as (x) can be 
deduced from p{x), using (x) = J_ xp{x) dx, and is zero as p{x) is symmetric. However 
we have seen that the positions of the minima in both Vcs{x) and V"s_ci(a;) are very close 
so the mean field values for {x)mf provide an estimate of the positions of the minima in 
Vcs{x). 

As (x) = in the exact solution for the symmetric model, a more interesting 
quantity to calculate is the root mean square deviation Ax, where (Aa;)^ = {{x — {x))"^). 
In Fig. [6] we give a plot of Ax deduced from p{x) for the model with U = and several 
values of a. 




U/ltA = 

©-OU/JtA = 2 



2 3 

a 



Figure 6. The root mean deviation Ax/xq as a function of a for the case U — 0, 
£f = (full curve), and U/nA — 2, Sf — ~U/2 (dashed curve). The dotted curve 
corresponds to ^/ottA/loq. 



It can be seen that Ax does vary significantly with a reflecting the fact that in the 
strong coupling regime the deviation Ax is determined by position of the minima in 
the double potential wells. As the mean field values (x)mf; were found to give a good 
estimate of these, and in the strong coupling regime (x)mf ~ ±A/a7rA/a;o, we compare 
in Fig. [6] this estimate with the calculated value of Ax. It is seen to provide a good 
estimate of Ax over the range of a shown. Also in Fig [6] we give the values of Ax for 
the symmetric model with U/ttA = 2. The tendency for the U term to suppress the 
fluctuations for the lower values of the coupling strength A, noted in Fig. [H is apparent 
but to have only a relatively minor effect in the stronger coupling range. 



2.2. Results for the asymmetric model 

With only a small degree of asymmetry, the form for p(x) changes quite dramatically in 
the strong interaction regime. This is because the doubly occupied impurity state is now 
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predominantly favoured. In Fig. [7]we show results for p{x) in a case with ef/nA = —0.5 
and f/ = for the same range of values of a as in Fig. [1] 
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Figure 7. The total probability distribution function p{x) for the oscillator 
displacement x in the ground state for a range of values of a with U = and 
ef/nA ^ -0.5. 



There is just a single narrow peak for p{x) in each case which shifts to slightly more 
negative values of x as a increases. In Fig. |8] we compare the semiclassical potential 
and Ves{x) as derived from equation ([7]) for the case a = 1. 




Figure 8. The effective potentials T4ff(a;) (full curve) as deduced from equation ([7]), 
compared with corresponding the semiclassical potentials V^_ci(2;) (dashed curve) for 
a = 2.0, Ef/nA = -0.5 and U = 0. 



Both potentials have an absolute minimum at x/xq = —1.29 though in the semiclassical 
case there is a secondary local minimum. This value of x/xq agrees with that predicted 
by equation (fl^ using the value derived from (fif) derived from the NRG calculation 
which gives x/xq = —1.2909. We can check the relation (IT^ by taking the average of 
X over the distribution p{x) and then compare with the result from equation (fT4|) using 



Probability distribution of fluctuations in Anderson- Holstein and Holstein- Hubbard Modelsll 

the NRG calculated (nj). The results for a range of values of a are shown as points 
(crosses and plus signs) in Fig. |9] (smaller values for f/ = 0, Ef/TiA = —0.5, and larger 
values U/nA = 2,ef = -U/2 - 0.05). 



cse-B-'" 




'■-© — o— — o 



— U ^ (mean field) 
■■■ U/ltA = 2 
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Figure 9. The values of {x)/xo for C/ = 0, ef/irA = -0.5, and U/nA = 2, 
Ef = —U/2 — 0.05, as calculated from the average of p{x) (crosses) and those deduced 
from equation (fT4|) (plus signs). The mean field results for C/ = correspond to the 
full curve. The corresponding results for the root mean deviation Ax/xq are also 
shown, for the case (7 = 0, Ef/irA — —0.5, (dashed curve and circles), and U/irA = 2, 
Ef = —U/2 (dotted curve and squares). 



For both sets of parameters, the results of the two calculations are in remarkably good 
agreement, giving the same values to at least five significant figures in all cases. The 
full curve in Fig. [9] corresponds to the mean field result (U = 0) for (x) /xq and can 
be seen also to be in good agreement with the exact results. Also shown in Fig. |9] is 
the root mean square Ax calculated from p(x). Though the average displacement {x) 
increases with increasing a it can be seen that Ax remains almost constant. In mean 
field theory Ax = 0, as in this approximation (x^) = (x)^. In the semiclassical approach 
one could estimate p(x) by solving the Schrodinger equation ([5]) with the potential 
V^_ci(a^) and using p(x) = |'?/'gs(x)p, where ipgs{x) in the ground state wavefunction, and 
then use the result to take an average of x^. We can, however, calculate it exactly 
in the limit of very weak and strong coupling limits. In the uncoupled case using 
the ground state wave function for the oscillator we find (Ax)^/xo = 1/2. In the 
strong coupling case with asymmetry we can take {fif) = 2 or (fif) = 0, and use a 
displaced oscillator transformation to new phonon creation and annihilation operators, 
a^, a^^-* = b^"^^ ± g/coo- The ground state |gs) then corresponds to the state a|gs) = 0, 
and in this state (b'^b) = g'^/ujQ, (x^) = Xq(1 + Ag'^/ujQ)/2. We have from equation IHM 
with {{uf — 1)) = ±1, (x)^ = 2g'^xl/u}Q, which gives again (Ax)^/xq = 1/2, so that we 
find Ax/xq = I/V2 in both limits. This agrees well with the results shown in Fig. [9] 
and Ax does not deviate much from this value over the whole range of a. As will be 
seen in the Sec. [3] for the case of the lattice model, the results will be different for Ax 
near the transition to a charge ordered state. 
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The corresponding values for the case with fZ/vrA = 2, ef/nA = —0.15 are shown 
in the same Fig. [9l The effect of finite U is to suppress the value of (x) for smaller 
values of a, but only has a marginal effect for a > 3 and has very little effect on Ax. 
Again there is five figure agreement in the the two calculations for (x); the one based 
on integrating over p{x) and the one using Eq. JHM . 

As the model in the limit Wq — t- cxd (A finite) corresponds to an Anderson model 
with an interaction term U — X. For U = and finite A, therefore, it becomes equivalent 
to a negative-f/ Anderson model. The symmetric model in the regime A/ttA ^ 1 has a 
Kondo effect due to charge rather than spin fluctuations. Introducing some asymmetry 
by changing Ef from the value for the symmetric case is equivalent to introducing a 
magnetic field in the Kondo case |11] . which for large fields suppresses the Kondo effect. 
We found that in using a value Ef/irA = —0.5 that the mean field estimate for (x) and 
the exact result were in good agreement. This is due to the fact that this degree of 
asymmetry corresponds to the Kondo case with a large magnetic field, and also because 
the value of uq used is much smaller than the bandwidth D. For a much smaller degree 
of asymmetry, which would correspond to a smaller 'magnetic' field, we should expect 
to find some limitations in the predictions from the semiclassical approximation. To 
examine this further we consider the case with a = 2 and Ef/iiA = —0.01. In Fig. [10] 
we show the semiclassical potential and Vef[{x) derived from ([7]). 
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Figure 10. The effective potentials T^h (x) (full curve) as deduced from equation ([7]), 
compared with corresponding the semiclassical potentials Vg-dix) (dashed curve) for 
a = 2.0, Ef/irA = -0.01 and [/ = 0. 



In this case we see that both potentials have two local minima. The absolute minimum 
in the two cases coincide at a value oi x/xq = —1.26, which corresponds to the mean field 
estimate of {x)/xo- The value obtained by averaging x over the distribution p{x), and 
from the formula (fl^ with the NRG value for (fif), both give {x)/xq = —1.09888. The 
average value in this case no longer coincides with the absolute minimum of the potential 
as the fluctuations to the local neighbouring minimum make a significant contribution. 
As the semiclassical equation ( ITT]) for (x) agrees with the exact one in equation ( TT^ . 
this difference arises from the fact that the semiclassical prediction for Uf disagrees with 
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the exact value of {rif)- The semiclassical prediction gives uj = 1.889 and the exact 
value from the NRG gives (hf) = 1.777, which explains the difference in the predictions. 
It is interesting to note, however, that the semiclassical prediction does coincide with 
the absolute minimum of the effective potential derived from the NRG results. 

The deviations from mean field theory become more marked for higher oscillator 
frequencies. In Fig. [11] we give a plot of {x)/xq for the case t/ = and Sf/ixA = —0.1 
taking the phonon frequency value Uq = 0.6. 
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Figure 11. The average oscillator displacement {x)/xo (full curve) as a function of a 
for the case [/ = 0, Ef/nA = —0.01, uiq ~ 0.6. The values calculated from the average 
over p{x) are indicated by crosses and those deduced from equation p4[) by plus signs. 
The dashed-dot curve gives the corresponding mean field results. Also shown are the 
results of the root mean square deviation Ax/xq (dashed curve and circles) . 



In this case, except for the small values of a, there is quite a discrepancy between 
the exact and mean field results. The corresponding results for the root mean square 
deviation Ax are also shown. The results for this quantity are very similar to those 
shown in Fig. |9l 

We can also examine the dependence of p{x) on the oscillator frequency Uq. In Fig. 
[T2] we show the change in form of p{x) as the frequency is increased for a fixed value of 
A in the strong coupling regime corresponding to a = 4.5 {U = 0, Ef = 0). 
We argued earlier that the minimum in the effective potential in the strong coupling 
regime occurs at a value of x ~ a/A/wq- We would expect the peak in p{x) to behave in 
a similar way, so for fixed A the peak positions should vary as I/cuq- This can be seen to 
be well satisfied in the results shown in Fig. [121 In the limit wq — )■ oo the double peak 
feature disappears entirely and p(x) becomes a delta function at x = 0. In this limit 
the mean field equation for Uf still has a broken symmetry solution for a > 0.5. This 
coincides with the static mean field solution for the Anderson model with U = —A, if one 
first performs a Hubbard- Stratonovich transformation to couple the auxiliary field x(r) 
solely to the impurity charge. It differs by a factor of 2 from the mean field theory of the 
Anderson model with U = —A, where the interaction term Und^^n^^i is approximated 
by U{nd,^{nd.x) +nd,|(nd,t) - (rid,t)('^d,;))- 
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Figure 12. The variation of p(x) with the frequency ojo in the strong coupUng case 
with a — 4.5. The value of xq is fixed and corresponds to l/y/UJo for ujo = 0.1. 



3. Local Fluctuations in the Holstein-Hubbard Model 

The states of broken symmetry predicted by the semiclassical/mean field theory for 
larger values of A cannot exist for the impurity Anderson- Holstein model; the symmetry 
has to be restored in the exact solution. Symmetry breaking, however, as a result of a 
phase transition can occur in a lattice model. To study p{x) in the neighbourhood of a 
phase transition we consider the Holstein-Hubbard model, described by the Hamiltonian, 

H = ^(ty4„Cj> + h.c.) + t/ ^ ^i^-fn-i^i (15) 



«J,o- 



+ UoJ2blbr + 9j2^h + bl) 



z^ 



Ui, 



1 



c|^ creates an electron at lattice site i with spin a, and b- a phonon with oscillator 
frequency loq, hi^„ = c\^Ci^fj. There is a coupling g to the local charge at each site, 
as in the Holstein model, and an on-site interaction U between spin-up and spin-down 
electrons, as in the Hubbard model. The hopping term tjj between orbitals localised on 
each site leads to a conduction band with a density of states Dffiuj) when g = U = 0. 
In the limit of infinite dimensions the model can be mapped into an effective Anderson- 
Holstein model, which can then be solved using the NRG. This has been described 
fully elsewhere [31 [2], and is known as the dynamical mean field theory (DMFT). 
For three dimensional systems, the mapping is only approximate. This approach is 
non-perturbative and can be applied in the strong interaction regime of these models 
to describe strong correlation effects, and indications are that it constitutes a good 
approximation when the self-energy of the electrons is local and a function of frequency 
only. 

There have been several applications of the DMFT method to study phase 
transitions in the Hubbard- Holstein model [T2l [T3| [HI [15]. There are various possible 
transitions to states of broken symmetry in this model; bipolaronic (BP), charge ordered 
(CO), antiferromagnetic (AFM) and the superconducting (SC) state. We restrict our 
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attention here to the case of half-filhng, so we do not include the superconducting case, 
which exists as a stable state away from half-filling. The transition first studied by the 
DMFT-NRG method for this model did not include the possibility of either charge order 
or antiferromagnetism |13 [ [l4j. There is, however, a metal- insulator transition from the 
normal state (N) to the bipolaronic state (BP), first studied for the model with [7 = 0, 
which occurs as the electron-phonon coupling A is increased at a critical value Ac- The 
transition also occurs in the model with f/ 7^ 0, at larger values of Ac, as the attractive 
term induced by A has to overcome the repulsion due to U{> 0). If the possibility of 
transitions to charge order (CO) and antiferromagnetism are included, then it has been 
found that antiferromagnetism occurs for U — X > and charge order for f/ — A < 
[161 [H]. 

For the DMFT-NRG calculations presented here we have taken a Bethe lattice form 
for the density of states Do{uj) of the conduction electrons, 

Doico) = ^V^^^^. (16) 

We choose a value t = 1 to set the energy scale in the following, which corresponds 
to a bandwidth W = At = 4. The physically relevant regime in the lattice case is for 
phonon frequencies small compared with the bandwidth. In most of the calculations in 
this section we take ojq = 0.6, which is small compared with the bandwidth of W but 
well away form the adiabatic limit. The distribution function for the local oscillator 
displacement p{x) was calculated from the DMFT-NRG density matrix using equation 
(ED. 



3.1. Model without long range order 

We consider first the results for the normal to bipolaronic transition, which are shown 
in Fig. [13] for U = 2 (left panel) and U = 5 (right panel) and different values of A. 





Figure 13. (Color online) The local probability distribution p{x) for U = 2 (right 
part) for the N state and U = 5 (left part) for various values of A for wq = 0.6. 
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In the left panel we can see for the case oi U = 2 how the probability distribution 
becomes broader as A is increased. As we do not allow for symmetry breaking here 
the system changes between zero occupation and double occupation with the associated 
oscillator fluctuations to minimise the energy. The situation is similar to the impurity 
case shown in Fig. [H and a two peak form develops in a similar way. However in this 
case when the two peak structure develops a gap also appears in the electron density of 
states, D{uj), signaling the transition to a insulating bipolaronic state. The correlation 
can be seen in the corresponding results for D{uj) shown in Fig. [Trover the transition 
regime. 




Figure 14. (Color online) The local electronic spectral functions I?(a;) in comparison 
for U = 2 (left panel) and U = 5 (right panel) for various A and wq = 0.6. 



This is in contrast to the impurity density of states shown in Fig. |3] where there is a 
shift of spectral weight from the region near u = to higher and lower values of u, but 
a narrow central peak at cj = remains. The narrow peak reflects the fact that there is 
no broken symmetry in the impurity case, and there are fluctuations between the two 
potential wells that restore the symmetry. In the right panel the results are shown for 
p{x) across the transition are shown for U = 5. There is a similar trend leading to a two 
peak structure in the bipolaronic phase, though at larger values of A due to the larger 
value of U. However, in this case there is an intermediate regime where p{x) has three 
maxima, which is very close to the metal-bipolaron transition, and the change is more 
marked, occurring over smaller range of A. 

Similar as in the impurity case, Eq. (JTj), we can compute the effective potential 
Ves{x), which is formed locally in the lattice model due to the electron-phonon coupling, 
from the probability distribution. This is shown in Fig. [15] for the same values as before, 
U = 2 (left panel) and U = 5 (right panel). 

In the case U = 2 on increasing A one sees that the potential becomes shallow and 
eventually develops two minima at finite ±Xm- Note that although there are two minima 
already for A = 2.8, fluctuations are sufficient to keep the system in a metallic state, as 
can be seen also from the spectral function in Fig. [Ml When the minima are deeper. 
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Figure 15. (Color online) The local effective potential V{x) ioi U ~ 2 (left part) for 
the N state for C/ = 5 (right part) for various values of A for ujq = 0.6. 



as for A = 3.2, the system is in the BP insulating state. The transition is continuous. 

In the case U = 5 the overall trend is similar, but larger values of A are required 
to induce the transition. Close to the transition we can find a structure of 3 local 
minima, where the one at x = is lifted upon increasing A. This is characteristic for a 
discontinuous transition which is expected to occur for larger values of U as discussed 
by KoUer et al. pH]. 

If we restrict to the pure Holstein model ([/ = 0), then it is also of interest here to 
study the quality of the semi-classical approximation. Similar as in the impurity case 
the potential can be calculated and one finds 



Vs-c\{x) 



-Do{mx)f-mxfDo{mx)) 



(17) 
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where ft{x) = fi — ^/2uJogx. The condition — °g^ = gives the mean field solutions and 
one can infer that at half filling, fi = 0, the potential has two minima if A > A™^ = nD/A. 
The value for this to occur in the DMFT with local quantum fluctuations is larger and 
also depends on Wq. For D = 2 one has, e.g., for Wq = 0.2 the value Ac ~ 1.75 and for 
cuo = 0.6 the value Ac — 2. For the first case, cuq = 0.2, we give a comparison of the 
effective potential obtained in DMFT calculations and the semiclassical approximation 
(fT7|) in Fig. [16] for two values of A. 

For the smaller value of A with the minimum at x = one finds a quite good agreement 
between the calculations. However, closer to the transition the results vary significantly. 
A = 1.6 is still a metallic solution with a narrow quasiparticle peak in DMFT. As 
A > A™^ the semiclassical approximation possesses two shallow minima for this case 
which can barely be resolved on the plot. The positions (1.16 semi-classical, 2.32 DMFT- 
NRG) differ significantly from the DMFT result, where fluctuations keep the state 
metallic. Better agreement between for the position of the minima in the semiclassical 
approximation and DMFT can again found in the insulating (BP) phase. 
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Figure 16. (Color online) Comparison of V{x) for DMFT-NRG calculation (full line) 
with the semi-classical approximation (dashed line) for ojq = 0.2. The values of ys_ci(a^) 
have been offset to make the comparison clearer. 



3.2. Model with long range order 

To consider charge ordered states, the lattice is divided into two sublattices denoted by 
A and B. Charge order develops when the occupation values {fii^A) 7^ {ni,B), and their 
difference divided by 2 can be taken as the order parameter $co- As {fii^A + ^j,b) = 2 at 
half filling, $co can be taken as l/2{{{hi^A) — !)• When charge order occurs p{x) shifts 
position to a displaced state appropriate to the local charge to minimise the energy. 
This can be seen in Fig. [T7] (left panel) where we show p{x) for U = 2 and values of A 
as charge order develops for A > f/. 




Figure 17. (Color online) The local probability distribution p{x) for [/ = 2 (left part) 
for the CO state for [/ = 5 (right part) for various values of A for wq = 0.6. 



In this case p{x) has a single peak which shifts and narrows slightly as A is increased. 
The same trend can be seen in the right panel of Fig. [T7]for the case U = 5, though the 
shifting and narrowing occurs more rapidly as A is increased. The shifting of a single 
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peak with increasing A is similar to the asymmetric impurity case shown in Fig. [TJ 
though the narrowing is an extra feature. 

The exact relation in equation (HM we derived earlier between the average 
displacement and the expectation value for occupation of the impurity site also holds for 
the lattice model. If we let {x) denote the value of (xi^A), then we have from equation 

m 









(18) 



so that the order parameter is directly proportional to (x). Again we can test this 
relation by calculating (x) from the average over p{x) and use the NRG results for $^05 
and find that it is satisfied very precisely. 

In Fig. [TS] we plot {x)/xq for various values of f/ as a function of A, which from 
equation ( TT8|) is proportional to the order parameter $co- 




Figure 18. (Color online) The expectation values 
function of A for uo — 0.6 in the CO state. 



for various values of [/ as a 



In the normal or antiferromagnetic state (x) =0, and the onset of the charge order can 
be seen clearly to occur for \ r^ U . The transition increases sharply with increasing f/, 
such that for [/ = 5 it is discontinuous. There is a similar trend in the impurity case 
shown in Fig. |6l but it is more marked in the lattice case. 

The mean square deviation (Ax)^ for the lattice coordinate, which can be deduced 
from the appropriate averages over p(x), is a measure of the fluctuations of the order 
parameter. In Fig. [19] (Aa;)^/xQ is plotted for the same set of parameters as for (x). 
This quantity is finite in the antiferromagnetic phase for X < U, and increases 
significantly as the transition is approached. The fluctuations can be seen to become 
much greater in the region of the transition for the intermediate values of U, most 
prominently near the point U = X = 4, where the transition changes from second to 
first order. This is a reflection of the fact that p{x) broadens at the transition and then 
narrows as A is increased further. Once the charge order has been well established the 
fluctuations then fall off rapidly to give {Ax^/xq = 1/2. Though the values of (Ax)^ 
away from the transition behave like the impurity case shown in Fig. [9l there is a very 



Probability distribution of fluctuations in Anderson- Holstein and Holstein- Hubbard Models20 



2 


-0-U=O 
-e-U=3 






1.5 


-i-U=4 
-^U=5 


i i 1 




(N O 
X 




' J\ i 




'^K 1 




\^ ff''\ /"' 




<l 


(y^^^^^^i: 


>K^i^ 


1 


0.5 








) 1 


2 3 4 5 



Figure 19. (Color online) The expectation values Ax'^ for various values of [/ as a 
function of A for ojq = 0.6 in the AFM state ior X < U and in the CO state for X > U. 



marked difference in the critical region especially for large U. An analysis of the effective 
potential is also possible for the CO case, but will be omitted here. For U = one can 
similarly derive a semiclassical potential. A detailed DMFT study of CO order in the 
Holstein model in the adiabatic limit, where also the probability distribution and the 
effective potential are discussed, has been given by Ciuchi et al. |18j. 



4. Conclusions 

We have shown how the reduced density matrix obtained from NRG calculations can 
provide physically relevant information about the local fluctuations. We have illustrated 
this by calculating the probability density function p{x) for the spatial coordinate x of 
the local oscillator in the impurity Anderson-Holstein and lattice Hubbard-Holstein 
models. This has enabled us to address a number of interesting questions. We have 
been able to see how the features in p{x) correlate with the features seen in the spectral 
density of the electrons as the interaction strength is increased. We have also deduced 



an effective potential Vcs{x), such that the wavefunction 



{X 



of the Schrodinger 



equation corresponds to p{x). This has enabled us to compare this potential with the 
one obtained from a semiclassical approximation, where the x-coordinate is treated as 
a classical variable, which is equivalent to the Born-Oppenheimer approximation. The 
results have provided a guide as to when the semiclassical approximation can be expected 
to give reliable results, and to clarify its limitations. 

We have also been able to compare the fluctuations of x in the impurity case with 
those in the lattice model in the various parameter regimes. For the normal state BP 
insulator transition we found a double well potential for weaker coupling and a structure 
with three local minima for stronger coupling. The semiclassical approach only gave 
a good description for weaker electron-phonon coupling. Allowing for the symmetry 
breaking in the lattice model, we have found that p{x) broadens in the critical region of 
the antiferromagnetic to charge order phase transition. The critical fluctuations become 
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particularly marked in the intermediate U regime near the point where the ground state 
transition changes from continuous to discontinuous behaviour. 

From a calculation of the reduced density matrix it is also possible to learn 
something about the local electronic fluctuations. If in equation (|2]) we integrate over 
the oscillator coordinate x, but do not carry out the sum over q and m^, then we have 
components p{q,mz) of the impurity reduced density matrix. From these, for example, 
we can deduce directly the impurity charge fluctuation, (An/)^ = {{fif — (n^))^), 

(An^)2 = 4p(2,0)+ Yl P(0,m,)-(2p(2,0)+ J^ p(0,m,))'.(19) 

m^=±l/2 m^=±l/2 

If the calculation of the reduced density matrix were to be terminated earlier, say at 
the neighbouring site to the impurity, local electronic fluctuations and near neighbour 
correlation functions could be deduced in a similar way. 
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